! VSemiletov, 08.2011
! Initialization for sendind all cubes to first processor which writes them to 'out' and 'outRes'
! Definition of:
! nCellsGlobal, nFacesGlobal,
! arrays for send and recv
subroutine OutputFromAllProcsInit()
use Data_cfoam
use Data_Mesh, only: nCells, nFaces, nCellsGlobal, nFacesGlobal
implicit none
include 'mpif.h'

    integer :: iProc

    if(nump.eq.1)then
        nCellsGlobal = nCells
        nFacesGlobal = nFaces
        return
    end if

    if(rank.eq.0)then
        allocate(recvCountsCell(nump),rDisplsCell(nump+1))
        allocate(recvCountsFace(nump),rDisplsFace(nump+1))
        recvCountsCell = 0
        recvCountsFace = 0
    else
        allocate(recvCountsCell(1),rDisplsCell(1))
        allocate(recvCountsFace(1),rDisplsFace(1))    
    end if

!    write(*,*) '@@@ MPI_Gather rank=', rank

    call MPI_Gather(nCells,1,MPI_INTEGER,recvCountsCell,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_Gather(nFaces,1,MPI_INTEGER,recvCountsFace,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

    if(rank.eq.0)then
        rDisplsCell(1) = 0
        rDisplsFace(1) = 0
        do iProc = 1,nump
            rDisplsCell(iProc+1) = rDisplsCell(iProc) + recvCountsCell(iProc)
            rDisplsFace(iProc+1) = rDisplsFace(iProc) + recvCountsFace(iProc)
!            write(*,*) 'Output from first proc: recvCountsCell(',iProc,')=', recvCountsCell(iProc), ' recvCountsFace=', recvCountsFace(iProc)
!            write(*,*) 'Output from first proc: rDisplsCell(',iProc,')=',    rDisplsCell(iProc+1),  ' rDisplsFace=',    rDisplsFace(iProc+1)
        end do
        
        nCellsGlobal = rDisplsCell(nump+1)
        nFacesGlobal = rDisplsFace(nump+1)
        
        write(*,*) 'nCellsGlobal=', nCellsGlobal, ' nFacesGlobal=', nFacesGlobal
        
        allocate(recvCell(nCellsGlobal))
        allocate(recvFace(nFacesGlobal))
    end if

end subroutine OutputFromAllProcsInit

! Send all cubes to first processor which writes them to 'out' and 'outRes'
subroutine OutputFromAllProcs()
use Data_cfoam
use Data_Mesh, only: nCells, nFaces, nCellsGlobal,nFacesGlobal
use DataTime, only: ttime, nstep, dt
implicit none
include 'mpif.h'

    integer :: i,j

    if(nump.eq.1)then
!        write(*,*) '1'
        write(32) ttime,0,nCells !,(U(1:3,i),p(i),T(i),i=1,nCells) !,E(i),rho(i)
        do j = 1,3
!            write(*,*) j+1
            write(32) (U(j,i),i=1,nCells) !,E(i),rho(i)
        end do
!        write(*,*) '5'
        write(32) (p(i),i=1,nCells) !,E(i),rho(i)
!        write(*,*) '6'
        write(32) (T(i),i=1,nCells) !,E(i),rho(i)
!        write(*,*) 'krutyak'
        
        write(outRestarts) ttime,nstep+1,dt,0,nFaces !,(Usnew(1:3,i),psnew(i),Rhosnew(i),i=1,nFaces) !,Esnew(i),Tsnew(i)
        do j = 1,3
            write(outRestarts) (Usnew(j,i),i=1,nFaces) !,Esnew(i),Tsnew(i)
        end do
        write(outRestarts) (psnew(i),i=1,nFaces) !,Esnew(i),Tsnew(i)
        write(outRestarts) (Rhosnew(i),i=1,nFaces) !,Esnew(i),Tsnew(i)
        write(outRestarts) 0,100 !,(utmp(1:8,i),i=1,100)
        write(outRestarts) (utmp(1:8,i),i=1,100)

        return
    end if

! Cells Output
    if(rank.eq.0)then
!        write(*,*) '@@@@ abracadabra=', ttime,0,nCellsGlobal
        write(outResults) ttime,0,nCellsGlobal ! 32
    end if
    do j = 1,3    ! Velocities Output
        do i = 1,nCells
!            U(j,i) = rank*100.0d0 + i
!            write(*,*) '@@@@ ', rank, ' U(',j,',',i,')=',U(j,i)
            VarMem(i) = U(j,i)
        end do
        
        call OutputArray(VarMem,nCells,recvCell,recvCountsCell,rDisplsCell,nCellsGlobal,outResults)
!        call MPI_Gatherv(VarMem,nCells,MPI_DOUBLE_PRECISION,recvCell,recvCountsCell,rDisplsCell, &
!                        MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
!        if(rank.eq.0)then
!!            do i = 1,nCellsGlobal
!!                write(*,*) '@@@@ Global U(',j,',',i,')=',recvCell(i)
!!            end do
!            write(outResults) (recvCell(i),i=1,nCellsGlobal) ! 32
!        end if
    end do
    ! Pressure Output
    call OutputArray(p,nCells,recvCell,recvCountsCell,rDisplsCell,nCellsGlobal,outResults)
!    call MPI_Gatherv(p,nCells,MPI_DOUBLE_PRECISION,recvCell,recvCountsCell,rDisplsCell, &
!                    MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
!    if(rank.eq.0)then
!        write(outResults) (recvCell(i),i=1,nCellsGlobal) ! 32
!    end if

    ! Temperature Output
    call OutputArray(T,nCells,recvCell,recvCountsCell,rDisplsCell,nCellsGlobal,outResults)
!    call MPI_Gatherv(T,nCells,MPI_DOUBLE_PRECISION,recvCell,recvCountsCell,rDisplsCell, &
!                    MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
!    if(rank.eq.0)then
!        write(outResults) (recvCell(i),i=1,nCellsGlobal) ! 32
!    end if
    
! Faces Output
    if(rank.eq.0)then
        write(outRestarts) ttime,nstep+1,dt,0,nFacesGlobal
!        write(*,*) '@@@@ abracadabra=', ttime,nstep+1,dt,0,nFacesGlobal
    end if
    do j = 1,3    ! Velocities Output
        do i = 1,nFaces
            VarMem(i) = Usnew(j,i)
        end do

        call OutputArray(VarMem,nFaces,recvFace,recvCountsFace,rDisplsFace,nFacesGlobal,outRestarts)
!        call MPI_Gatherv(VarMem,nFaces,MPI_DOUBLE_PRECISION,recvFace,recvCountsFace,rDisplsFace, &
!                        MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
!        if(rank.eq.0)then
!            write(outRestarts) (recvFace(i),i=1,nFacesGlobal) !,E(i),rho(i)
!        end if
    end do
    ! Pressure Output
    call OutputArray(psnew,nFaces,recvFace,recvCountsFace,rDisplsFace,nFacesGlobal,outRestarts)
!    call MPI_Gatherv(psnew,nFaces,MPI_DOUBLE_PRECISION,recvFace,recvCountsFace,rDisplsFace, &
!                    MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
!    if(rank.eq.0)then
!        write(outRestarts) (recvFace(i),i=1,nFacesGlobal) !,E(i),rho(i)
!    end if
    
    ! Density Output
    call OutputArray(Rhosnew,nFaces,recvFace,recvCountsFace,rDisplsFace,nFacesGlobal,outRestarts)
!    call MPI_Gatherv(Rhosnew,nFaces,MPI_DOUBLE_PRECISION,recvFace,recvCountsFace,rDisplsFace, &
!                    MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
    if(rank.eq.0)then
!        write(outRestarts) (recvFace(i),i=1,nFacesGlobal) !,E(i),rho(i)
        write(outRestarts) 0,100 !,(utmp(1:8,i),i=1,100)
        write(outRestarts) (utmp(1:8,i),i=1,100)
    end if

end subroutine OutputFromAllProcs

subroutine OutputArray(Arr,nElems,recv,recvCounts,rDispls,nElemsGlobal,out)
use Data_cfoam, only: nump, rank, ierr
implicit none
include 'mpif.h'

    integer,intent(in) :: nElems,nElemsGlobal, out
    integer,intent(in) :: recvCounts(nump),rDispls(nump+1)
    real(8),intent(in) :: Arr(nElems)
    real(8),intent(out) :: recv(nElemsGlobal)

    integer :: i
    
    recv = 0.0d0

    call MPI_Gatherv(Arr,nElems,MPI_DOUBLE_PRECISION,recv,recvCounts,rDispls, &
                    MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
    call MPI_Barrier(MPI_COMM_WORLD,ierr)

    if(rank.eq.0)then
        write(out) (recv(i),i=1,nElemsGlobal)
    end if
    
end subroutine OutputArray

subroutine WriteMeshCellFaceCenters()
use Data_cfoam, only: recvCell,recvCountsCell,rDisplsCell, &
                      recvFace,recvCountsFace,rDisplsFace, &
                      rank, nump, ierr
use Data_Mesh, only: nCells, nFaces, meshCellCenters, meshFaceCenters, nCellsGlobal, nFacesGlobal
implicit none

    real(8),allocatable :: CellsMem(:), FacesMem(:)

    integer :: i,j

    if(rank.eq.0)then
        write(*,*) 'meshCellCenters, meshFaceCenters are being written ...'

        open(11,file='outMeshCellCenters',form='unformatted',iostat=ierr)
        if(ierr.ne.0) then
   	        stop '@@ Fatal Error: file outMeshCellCenters has not been opened!'
   	    end if
        open(12,file='outMeshFaceCenters',form='unformatted',iostat=ierr)
        if(ierr.ne.0) then
   	        stop '@@ Fatal Error: file outMeshFaceCenters has not been opened!'
   	    end if
   	end if

    if(nump.eq.1)then
        write(11) 0d0,3,nCells
        do j = 1,3
            write(11) (meshCellCenters(j,i),i=1,nCells)
        end do
        write(12) 0d0,3,nFaces
        do j = 1,3
            write(12) (meshFaceCenters(j,i),i=1,nFaces) !,E(i),rho(i)
        end do
    else
        allocate(CellsMem(nCells))

        if(rank.eq.0)then
            write(11) 0d0,3,nCellsGlobal
        end if

        do j = 1,3
            do i = 1,nCells
                CellsMem(i) = meshCellCenters(j,i)
            end do

            call OutputArray(CellsMem,nCells,recvCell,recvCountsCell,rDisplsCell,nCellsGlobal,11)
        end do
        
        deallocate(CellsMem)

        ! meshFaceCenters output
        allocate(FacesMem(nFaces))

        if(rank.eq.0)then
            write(12) 0d0,3,nFacesGlobal
        end if
        do j = 1,3
            do i = 1,nFaces
                FacesMem(i) = meshFaceCenters(j,i)
            end do

            call OutputArray(FacesMem,nFaces,recvFace,recvCountsFace,rDisplsFace,nFacesGlobal,12)
        end do

        deallocate(FacesMem)
    end if

    close(11)
    close(12)

    if(rank.eq.0)then
        write(*,*) 'meshCellCenters, meshFaceCenters have been written successfully!'
    end if

end subroutine WriteMeshCellFaceCenters
